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On the thermo-elastostatics of heterogeneous materials 

II. Analyze and generalization of some basic hypotheses and propositions 



Abstract One considers linearly thermoelastic composite media, which consist of a homogeneous matrix con- 
taining a statistically homogeneous random set of ellipsoidal uncoated or coated inclusions. Effective properties 
(such as compliance and thermal expansion) as well as the first statistical moments of stresses in the phases are 
estimated for the general case of nonhomogeneity of the thermoelastic inclusion properties. At first, one shortly re- 
produces both the basic assumptions and propositions of micromechanics used in most popular methods, namely: 
effective field hypothesis, quasi-crystallite approximation, and the hypothesis of "ellipsoidal symmetry". The ex- 
plicit new representations of the effective thermoelastic properties and stress concentration factor are expressed 
through some building blocks described by numerical solutions for both the one and two inclusions inside the 
infinite medium subjected to both the homogeneous and inhomogeneous remote loading. The method uses as a 
background the new general integral equation proposed in the accompanied paper and makes it possible to aban- 
don the basic concepts of micromechanics mentioned above. The results of this abandonment are quantitatively 
estimated for some modeled composite reinforced by aligned continously inhomogeneous fibers. Some new effects 
are detected that are impossible in the framework of a classical background of micromechanics. 
Keywords: A. microstructures, B. inhomogeneous material, B. elastic material. 



1. Introduction 

The prediction of the behavior of composite materials in terms of the mechanical properties of 
constituents and their microstructure is a central problem of micromechanics, which is evidently reduced 
to the estimation of stress fields in the constituents. Appropriate, but by no means exhaustive, references 
for the estimation of effective elastic moduli of statistically homogeneous media are provided by the 
reviews [1-7]. It appears today that variants of the effective medium method by Kroner [8] and by Hill 
[9], and the mean field method [10], [11] are the most popular and widely used methods. Recently a new 
method has become known, namely the multiparticle effective field method (MEFM) was put forward 
and developed by the author (see for references Buryachcnko [6]). The MEFM is based on the theory of 
functions of random variables and Green's functions. Within this method one constructs a hierarchy of 
statistical moment equations for conditional averages of the stresses in the inclusions. The hierarchy is 
then cut by introducing the notion of an effective field. This way the interaction of different inclusions is 
taken into account. Thus, the MEFM does not make use of a number of hypotheses which form the basis 
of the traditional one-particle methods. 

However, a diversity of micromechanical methods and their specific formulations astonish our imagi- 
nation only at first glance. We will see that most popular methods are based just on a few basic concepts 
of micromechanics. Effective field hypothesis is apparently the most fundamental, most prospective, and 
most exploited concept of micromechanics. This concept has directed a development of micromechanics 
over the last sixty years and made a contribution to their progress incompatible with any another con- 
cept. The notion of an effective field in which each particle is located is a basic concept of such powerful 
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methods in micromechanics as the methods of self-consistent fields and effective fields (for references see 
[6], [7], [12]. The idea of this concept dates back to Mossotti [13] and Clausius (in the dielectric context), 
Lorenz (in the refractivity context), and Maxwell (in the conductivity context). Markov [14] and Scaife 
[15] presented comprehensive reviews of the 150-year history of this concept accompanied by some famous 
formulae with extensive references. Mossotti [13] (especially Clausius) pioneered the introduction of the 
effective field concept as a local homogeneous field acting on the inclusions and differing from the applied 
macroscopic one. Among a few hypotheses used by Mossotti [13], the most important one was in fact 
the quasi-crystalline approximation proposed 100 years later by Lax [16] in a modern concise form. The 
concept of the effective field in combination with subsequent assumptions was introduced in a modern 
formalized form in the physics of multiple scattering of waves (see, e.g., [16-18]). Walpolc [19] pioneered 
the application of the concept to the static of composites under the name uniform image field. Effective 
field technique was intensively applied in micromechanics of random and periodic structure composites 
(for references see, e.g., [6], [7]) as well in micromechanics of multiple interacting cracks under the name 
traction or pseudo-load [20]. Buryachenko and Rammerstorfer [21] has drawn the conclusion that the 
effective field concept is used (either explictly or implicitly) in most popular methods of micromechan- 
ics such as, e.g., the effective medium method and their modifications, differential scheme, Mori-Tanaka 
method, and, needless to say, the MEFM. 

The idea of effective field and quasi-crystalline approximation were added by the hypothesis of 
"ellipsoidal symmetry" for the distribution of inclusions attributed to Willis [22]. As a tool for concrete 
applications of the concepts mentioned, the Eshelby [23] solution was used although the Eshelby's theorem 
has a fundamental conceptual sense (it will be shown in the current paper) rather than only an analytical 
solution of some particular problem for the ellipsoidal homogeneous inclusion. All these concepts creating 
the framework and background of modern statistical analytical micromechanics were transformed by 
the use of both the additional assumptions and sophisticated analytical and numerical tools to a few 
particular methods. However, we will show in this paper that the effective field hypothesis (also called 
the hypothesis Hla is a central one and other concepts play a satellite role providing the conditions for 
application of the effective field hypothesis. Moreover, we will show that all mentioned hypotheses are 
not really necessary and can be relaxed. 

The outline of the study is as follow. In Section 2 we recall the basic concepts defining the background 
of micromechanics. The interconnection between the different concepts and their essence are established. 
In Section 3 the auxiliary problem for one inclusion in the infinite matrix is presented for a general 
remote loading. The new general integral equation obtained in an accompanying paper by Buryachenko 
[24], henceforth referred to as (I), is presented in Section 4 through the operator forms of the particular 
solutions for both one and two interacting inclusions. This equation is solved by the iteration method in 
the framework of the quasi-crystallite approximation but without basic hypotheses of classical microme- 
chanics such as both the effective field hypothesis and "ellipsoidal symmetry" assumption. In Section 5 we 
qualitatively explain the advantages of the new approach with respect to the classic ones and demonstrate 
the corrections of popular propositions obtained in the framework of the old background of micromechan- 
ics. Quantitative estimations of results of the abandonment of the central hypothesis Hla are presented 
in Section 6. 

2. Preliminaries. Basic assumptions and propositions of micromechanics 

2.1 General integral representations and notations 

For the sake of brevity of the current presentation, the basic equations of thcrmoelasticity, the homo- 
geneous boundary conditions (2.51), statistical description of the composite microstructure, assumption, 
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and notations exploited in the current paper are presented in the accompanied paper by Buryachcnko 
[24] and the interested reader is referred to this publication, henceforth referred to as (I). 

In this section we will shortly reproduce both the basic assumptions and propositions of micromc- 
chanics in the form adopted for subsequent presentation. In most detail we will consider the mentioned 
concepts as applied to the MEFM based on some mathematical approximations for solving the infinite 
systems of integral equations involved, although other methods exploiting these concepts will also be 
discussed. 

For simplicity, we will consider only statistically homogeneous media (described, as a particular case, 
in Section 2.2 in I) subjected to the homogeneous boundary conditions (2.51). If clastic properties of the 
comparison medium and matrix coincide (3.201) then the known general integral equation in terms of 
stresses (see, e.g., [6], I) 

o-M = (o-)+y[r( X -y)77(y)-(r(x-y)77)(y)]dy, (2.1) 

t(x) = (a)+jT{x-y)[ V (y)-( V )(y)]dy, (2.2) 

can be much easier to solve because the stress-strain fields can be studied inside the inhomogeneities but 
not in the matrix; here rj = M.icr+f3 1 and T are the strain polarization tensor (3.141) and the Green stress 
tensor (3.51), respectively. Buryachcnko [6, I] proved that for no long-range order assumed, and for x G w 
considered in Eqs. (2.1) and (2.2) and removed far enough from the boundary r (a <C |x — y|, Vy G T), 
the right-hand side absolutly convergent integrals in (2.1) and (2.2) do not depend on the shape and 
size of the domain w, and they can be replaced by the integrals over the whole space R d . With this 
assumption we hereafter omit explicitly denoting R d as the integration domain in the equations. The 
new exact equation (2.1) forming a new background of micromechanics yields the known approximate 
one (2.2) only at some additional assumptions [(3.231) or (3.241), see for details I]. 

The solution of Eqs. (2.1) and (2.2) by the use of different assumptions provides the estimations of 
both the effective compliance M* and the effective eigenstrains (3* governed by the overall constitutive 
relation 

<e)=M»+/3* (2.3) 

and defined by general relations 

M* = M (0) + (MiB*), (2.4) 

/3*=^°) + (B* T /3 1 ), (2.5) 

where B* = B*(x) (x G v) is a local stress concentration tensor in the inhomogeneities obtained under 
pure mechanical loading (/3 = 0) 

er(x) = B*(x)(er) for x G v. (2.6) 



Analysis of Eqs. (2.1)-(2.6) leads to the universally accepted 

Proposition 1. For statistically homogeneous media subjected to the homogeneous boundary conditions, 
linear elastic effective properties M* and f3* depend only on stress distributions inside the inhomogeneities 
v but not inside the matrix V^ ' . 

Let the inclusions v\, . . . , v n be fixed and we define two sorts of effective fields <Ti(x) and <Ti....,n(x) (i = 
1, . . . , n; x G vi, . . . , v n ) by the use of the rearrangement of Eq. (2.1) [or (2.2)] in the following form (see 
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for the earliest references of related manipulations [6]): 
<t(x) = <fi(x) + J r(x - y)y i (y)r ? (y)dy, 
ff»(x) = CTi,...,„(x)+^ A r(x-y)^(y)77(y)dy ! 

2ri,.„,„(x) = (<T)(x)+ / {r(x-y)77(y)y(y|;u 1 ,x 1 ;...;t;„,x„)-(r(x-y)T ? )(y)}dy, (2.7) 

for x G Uj, i = 1,2, here V(y|;«i,Xi; . . . ;v n ,x n ) is a random indicator function of inclusions 

x G v under the condition that Xj ^ Xj if i ^= j (i,j = 1, . . . ,n). Then, considering some conditional 
statistical averages of the general integral equation (2.1) leads to an infinite system of new integral 
equations (n = 1,2,...) 

(cr|w 1 ,xi;...;u„,x„)(x) - ^ / T(x - y) (K t (y)j7|vi, x x ; . . . ; u„, x„)(y)dy 

i= 1 

= (<x)(x) + J {r(x-y)(r 7 |;« 1 ,x 1 ;...; U „,x„)(y)-(r(x-y)r 7 )(y)}dy. (2.8) 

Since x G «i, . . . , v„ in the n-th line of the system can take the values of the inclusions Vi, . . . , v n , the 
n-th line actually contains n equations. The definitions of the effective fields CTj(x), 0"i,2,....n( x ) a s well as 
their statistical averages (<Fj)(x), (cti 2,... n)( x ) are nothing more than notation convenience for different 
terms of the infinite systems (2.7) and (2.8), respectively. The physical meaning of these fields and their 
graphic illustrations are presented in Ref. [6]. 

2.2 Approximate effective field hypothesis 

In order to simplify the exact system (2.8) we now apply the so-called effective field hypothesis which 
is the main approximate hypothesis of many micromechanical methods: 

Hypothesis la, Hla). Each inclusion Vi has an ellipsoidal form and is located in the field (2.I2) 

Wi(y) = W{xi) (y e ih) (2.9) 

which is homogeneous over the inclusion Vi. 

In some methods (such as, e.g., the MEFM) this basic hypothesis Hla is complimented by a satellite 
hypothesis: 

Hypothesis lb, Hlb) The perturbation introduced by the inclusion V\ at the point y ^ v $ is defined by 
the relation 

J T(y - x)K(x)rj(x)dx = ViTifr - xjrii. (2.10) 

Hereafter ■q i = (rj (x) Vi (x) ) (j) is an average over the volume of the inclusion Vi (but not over the ensemble) , 
((•)>* = (((•))(<))) and ( x e v h y G Vj) 



(2.11) 



where the tensor is associated with the well-known Eshelby tensor by Sj = I M^Qj. For a 
homogeneous ellipsoidal inclusion vi the standard assumption (2.9) (see, e.g., [6], [7]) yields the assumption 
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(2.10), otherwise the formula (2.10) defines an additional assumption. The tensors Ty (x^ — x^-), proposed 
by Willis and Acton [25] for identical spherical inclusions have an analytical representation for spherical 
inclusions of different size in an isotropic matrix (see for references [6] ) regardless of whether the inclusions 
are coated or uncoated. 

According to hypothesis Hla and in view of the linearity of the problem there exist constant fourth 
and second-rank tensors B.j(x), R 2 (x) and Cj(x), Fj(x), such that 

<x(x) = B 4 (x)<r(x ? ) + Ci(x), ViT](x) = Ri(x)w(xi) + Fi(x), x e v i} (2.12) 

where v l C v w and R 4 (x) = w l M^ ) (x)B,(x), F;(x) = w 4 [M^ i) (x)C, i (x) +/3 x (x)]. According to Eshelby's 
[23] theorem there are the following relations between the averaged tensors (2.12) = WiQ~ 1 (I — Bj), 
Fj = — TJjQ J ~ 1 Cj, where gj = (g(x))/j) (g stands for B,C,R, F). It should be mentioned that the field 
CT(Xj) can vary with the location of the center Xj of the inclusion considered, but the field W(y) (y G u;) 
is homogeneous over the inclusion Wj. Because of this the application of Eshelby's theorem is correct. 

For example, for the homogeneous ellipsoidal domain m with 

M^(x) = = const, pf{x) = Pi* = c °nst at x G v h (2.13) 

we obtain 

B, = (i + Q.M^)" 1 , C< = -BiQi/3®. (2.14) 

In the general case of coated inclusions Vi, the tensors Bj(x) and Cj(x) can be found by the transformation 
method by Dvorak and Benveniste [26] (see for references and details [6], [27]). 

Using hypothesis HI (combining the hypotheses Hla and Hlb), the system (2.72) for k fixed 
inclusions with fixed values CTi,...,fc(x) (x G Uj, i = l....,fc) on the right-hand side of the equations 
becomes algebraic when the solution (2.12) for one inclusion in the field <f(Xj) (i = 1, . . . , k) is applied 

k 

R.Tri x, : +Fi = J2 Z,.,{k,,T; ,.;x,: • K,}. (2.15) 

i=i 

where the matrix Z~ x has the elements (Z _1 )y 

(Z-% = 1% - (1 - ^/iR/r.Jx, - x,), ■: /../ 1 n). (2.16) 

2.3 Closing effective field hypothesis and effective properties 

Different methods can be employed (see for details [6]) to truncate the hierarchy (2.8) considered as a 
system of coupled equations. One begins with the last hierarchy item which has the most heterogeneities 
held fixed, because this equation does not depend on the other. The solution obtained presents the forcing 
term in the next equation up the hierarchy. The unconditionally average field is finally obtained by going 
step by step up the hierarchy. For termination of the hierarchy of statistical moment equations (2.8) we 
will use the closing effective field hypothesis: 

Hypothesis 2a, H2a) For a sufficiently large n, the system (2.8) is closed by the assumption 

(<7i j n+i(x))j = (5"i i ... ! n(x))i, where the right-hand-side of the equality does not contain the index 

j ^ i (i = 1, . . . , n; j = 1, . . . , n + 1; x g Vi). 
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The hypothesis H2a rewritten in terms of stresses <r(x), (x S Vi) is a standard closing assumption 
(see e.g. [1], [28]) degenerating to the "quasi-crystalline" approximation [16] at n = 1 (see for analysis 
also Subsection 2.4). 

In the framework of the hypothesis HI, substitution of the solution (2.12), and (2.15) (at k = 2) for 
binary interacting inclusions into the first equation of the system (2.8) at n = 1 and at the effective field 
hypothesis H2a) with the first order approximation: 



(<Tj jg (x))j = (ct(x))j = const, x e vj(j =i,q). 



(2.17) 



leads to the solution (x € Vi] 



(<r) i (x) = y > T 9 (x-xJZ ?l <^(i; 9 ,xJ;v l ,xOdx 9 (R i (er} l +F ; ) 

T q {x-x q )Z qq tp(v q ,x q \;v i ;x i )-T{x-x q )n^] (R q (a) q + F q )dx q , (2.18) 



where the matrix elements Z q i, Z qq are nondiagonal elements and diagonal ones of the binary inter- 
action matrix Z (2.16) for the two inclusions v q and Vi\ hereafter the conditional probability density 
(p(v q , x,j|; Vi] Xi) and probability density ip(v q ,x q ) = are described in (I). Averaging the result ob- 
tained (2.18) over the inclusion Vi yields the final representation for both the statistical average stress 
field and effective properties 



iY 



( < r) i (x) =B l (x)R- 1 [^Y. y (R J ( ( r) 



Ci(x), 



N 



N 



M* = M (0) + J2 Y,,R ; // ' . [3* = /3 (0) + Y v F J 



(2.19) 
(2.20) 



where the matrix Y determines the action of the surrounding inclusions on the considered one and has 
an inverse matrix Y _1 given by 



(Y )ij — 5ij 

— R; 



JV 

I-R 4 ^ / T ig (xi - x q )Z qi ip(v q ,x. q \;v i ,-K i )dx. (i 

q=l J 



T ig (x,; - X. q )Z qq (p(v q ,-X q \;Vi,-X4) -Tj(Xi - X,)' 



dx q . 



(2.21) 



Buryachenko [6] demonstrated that the MEFM includes in particular cases the well-known methods of 
mechanics of strongly heterogeneous media (such as the effective medium and the mean field methods). 

2.4 Quasi-crystalline approximation 

Hypothesis H2a rewritten in terms of stresses <r(x) degenerates to the "quasi-crystalline" approxi- 
mation by Lax [16] which in our notations has two equivalent forms 

Hypothesis 2b, H2b, "quasi-crystalline" approximation. It is supposed that the mean value of 
the effective field at a point x€Oi does not depend on the stress field inside surrounding heterogeneities 

Vj ^ Vf. 

(^(x^UjjXjju^Xj) = (Wi), xeu,, or Zij=lSij. (2.22) 
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Therefore, the matrix Y 1 can be reduced to (see Ref. [29]) 

(Y-%- = 1% - R. t J fr i3 -(xj - XjMvj^Vi,*) - Ti(xj - x^n^dxj. (2.23) 

The principal difference between the hypotheses (2.17) and (2.22) are discussed in Chapter 9 in [6]). 

Note. It should be mentioned that the hypotheses H2a and H2b are not conceptually dependent on the 
hypothesis HI and can be applied in general case even if the hypothesis HI is violated (see for details 
Subsection 5.3). 

2.5 Hypothesis of "ellipsoidal symmetry" of composite structure 

To make further progress, the hypothesis of "ellipsoidal symmetry" for the distribution of inclusions 
attributed to Willis [22] (see also Khoroshun [28], [31], Buryachenko and Parton [29], Ponte Castaneda 
and Willis [32]) is widely used: 

Hypothesis 3, H3, "ellipsoidal symmetry". The conditional probability density function <p(vj,Xj 
; Uj, Xi) depends on Xj — Xi only through the combination p = |(aiy) (xj — Xj)|: 

<f(vj,x.j \]Vi,Xi) = h(p), p =| (a?-) _1 (xj - x,) | (2.24) 

where the matrix (aq^) (which is symmetric in the indexes i and j, a?- = a° iy ) defines the ellipsoid 
excluded volume = {x : |(a^) _1 x| 2 < 1}. 

A pair distribution function has "ellipsoidal symmetry" but with an ellipsoid shape differing from the 
one that defines the inclusion shape. Although the assumed statistics may not be exactly realized in any 
particular composite, the results of effective moduli estimations are explicit and simple to use. It is crucial 
for the analyst to be aware of their reasonable choice of the shape of "ellipsoidal" spatial correlation of 
inclusion location (see Chapter 18 in Ref. [6]). For spherical inclusions the relation (2.24) is realized for a 
statistical isotropy of the composite structure. It is reasonable to assume that (a^-) identifies a matrix 
of affine transformation that transfers the ellipsoid vfj being the "excluded volume" ( "correlation hole" ) 
into a unit sphere and, therefore, the representation of the matrix Yjj can be simplified: 

(¥"%•= Mij-RiQ?., (2.25) 

where Q°- = Q(«y) is a constant for the ellipsoidal domain with the indicator function V®j. For the 
sake of simplicity of the subsequent calculation we will usually assume that the shape of "correlation 
hole" u° does not depend on the inclusion vf. u° = and Q° = Q" = Q(w°). 

Substitution of the representation (2.25) into Eqs. (2.19) and (2.20) completes the problem of effective 
properties estimations. The hypothesis H3 is widely used for micromechanical structures described also 
by the indicator function V(x|; v\, Xi; ...;«„, x n ) (see, e.g. Willis [1], [22]; Ponte Castaneda and Willis 
[32]) rather than only by the conditional probability density (p(v q ,x q |;i>j,Xj) (2.24). 

Note. A popular point of view is that the hypothesis of "ellipsoidal symmetry (2.24) is exploited just 
for some simplification of the representation (2.23) reduced to (2.25). However, we will demonstrate in 
Section 5 that the destination of the hypothesis H3 is more fundamental and directed towards providing 
of conditions for applying of the hypothesis HI. The use of the satellite hypothesis H3 has no sense 
without the hypothesis HI. 

Proposition 2. If the hypotheses HI, H2b, H3 hold for the statistically homogeneous medium and 
homogeneous boundary conditions then the effective properties M* and f3* do not depend on the size of 
the correlation hole v® and the conditional probability density (2.24). 
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We can reach this conclusion by simple analyses of the final representations (2.19), (2.20), (2.23), 
and (2.25) as well as of analogous representations obtained by other methods in the framework of the 
hypotheses HI, H2b, and H3 (see, e.g., [1], [22], [28]). Here we keep in mind, first of all, the MEF which is 
equivalent for aligned identical ellipsoidal inclusions to the Mori-Tanaka [10] method. Moreover, Markov 
[33] demonstrated that for homogeneous ellipsoidal inclusions, the estimations by the MEF coincide with 
the variational estimates obtained by Ponte Castaneda and Willis [32] who also exploited the hypotheses 
HI, H2b, and H3. It will be demonstrated in Section 6 that satisfiability of the hypotheses H2b and 
H3 without HI can lead to dependence of M* and (3* on both the size of the correlation hole and 
the binary correlation function. 

3. A single inclusion subjected to inhomogeneous prescribed effective field 

3.1 Operator representation for solution obtained by the VIE method 

In the current section we will present a slightly modified solution of a satellite problem (sec [6] where 
additional references can be found) which is adapted for estimation of effective properties of composites 
in Section 6. Namely, let the inclusions Vi be fixed and loaded by the inhomogeneous effective field <x.;(x): 



<r(x) = <Ti(x) + J r(x - y)Vi(y)ri(y)dy, (3.1) 



A known quadrature method for obtaining an approximate solution of Eq. (3.1) is to evaluate the volume 
integrals with a Gauss quadrature formula. Then the corresponding equations at the Gauss points will 
contain a singular term that results when the field point x and source point y coincide: x = y (3.1). The 
difficulties with the troublesome singularities can be avoided if a rearrangement of Eq. (3.1) is performed 
in the spirit of a subtraction technique used in the modified quadrature method (see, e.g., [34]) 

tj(x) - M 1 (x)ct(x)+/3 1 (x)+Mi(x) /V i0 (y)r(x-y)d y77 (x) 

+ Mx(x) /'r(x-y)[T7(y)-y i o(y)T?(x)]dy, x G u,;, (3.2) 

where Vi, perhaps, is not an ellipsoid. The equation is valid for any domain x G i>io with the indicator 
function V^o(x). We assume that the first integral in (3.2) is easily computable 

- Q i0 (x) = J K i0 (y)r(x - y)dy. (3.3) 

In a general case of the inclusion shape Vi the function Q i0 (x) can be found numerically, e.g. by finite 
element analysis (see, e.g., [35]). For ellipsoidal domain vto D Vi, the first integral on the right-hand-side 
of (3.2) is known and is associated with the well-known Eshelby tensor by (x G Vi C vto, y G Vio) 

S i0 = I - M (0) Q ?;o , Q i0 = -«io<r(x - y))° = const. (3.4) 

Hereafter 

Si = (g(y))i = vr 1 J g(y)^(y)dy, g ?0 = (g(y)>° = (^o)" 1 / g(y)Ko(y)<fy (3.5) 



denotes averaging of some tensor g(y) over the volume of the regions y e u, C vm and y G v^, respec- 
tively; in so doing, x G vi is fixed. The assumption of an ellipsoidal shape of the domain Vio was used 
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only to obtain analytical representation of the integral (3.4). This is because the tensor (r(x — y))io is 
homogeneous for xGv;, y G Vio for an ellipsoid. For nonellipsoidal inclusions v n ~ e one could assume that 
in some parts of the region C Vi the properties Mi(x) = 0, /3i(x) = 0, i.e. it is sufficient to replace a 
real nonellipsoidal inclusion v n ~ c = Vi\vf by a fictitious ellipsoid (with smallest possible volume) and call 
it the inclusion m with a "coating" v^. In so doing, at the estimation of the second integral in (3.2) we 
keep in mind that / r(x — y)(V;(y) — Vio(y))dy = Q;o — Q;( x ) (xBj), With the nonessential restriction 
on the shape of the inclusion Vi mentioned above, we can consider without loss of generality an ellipsoidal 
inclusion = Vi] in so doing Mi(y) = 0, /?i(y) = at y G u| C Vi. Then Eq. (3.2) can be rewritten in 
the equivalent compact form 

V(x) = »7;W + y K ( (x,y)[r7(y) -r)(x)]dy, xBj, (3.6) 

where Vi( x ) = Ej(x)CT(x)+Hj(x), (x e w,) is called the effective strain polarization tensor in the inclusion 
Vi, and (no sum on i) 

Ki(x,y) = E l (x)r(x-y)^(y), (3.7) 
E 4 (x) = M 1 (x)[I + Q i0 (x)M 1 (x)]- 1 , (3.8) 
Hi(x) = [I + M 1 (x)Q 40 (x)]- 1 /3 1 (x). (3.9) 



We rewrite Eq. (3.6) in symbolic form 
where 



r 7 =^ + /C i » 7 , (3.10) 

(K i r,)(x) = J 7Q(x,y),7(y)dy (3.11) 
defines the integral operator 7C; with the kernel formally represented as 

7Q(x, y) - K s (x, y) - «5(x - y) J ^(z)K. i (x, z)dz. (3.12) 

We formally write the solution of Eq. (3.10) as 

xi = Ci*rj it (3.13) 

where the inverse operator C{ = (I — K,^ 1 will be constructed by the iteration method based on the 
recursion formula 

V [k+1] = % + KiV lk] (3.14) 

to construct a sequence of functions {j7^} that can be treated as an approximation of the solution of 
Eq. (3.10). We presented the point Jacobi (called also Richardson and point total-step) iterative scheme 
for ease of calculations. The details of the real iteration method used for the solution of Eq. (3.14) will 
be presented in Section 6. Usually the driving term of this equation is used as an initial approximation: 

r,[°](x)=W i (x), (3.15) 

which is exact for a homogeneous ellipsoidal inclusion subjected to remote homogeneous stress field 
<x(x) = W = const. The sequence {17^} (3.14) with arbitrary continuous rj^(x) converges to a unique 
solution 77 if the norm of the integral operator 7Q turns out to be small "enough" (less than 1), and 
the problem is reduced to the computation of the integrals involved, the density of which is given. In 
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effect the iteration method (3.14) transforms the integral equation problem (3.14) into the linear algebra 
problem in any case. 

We will introduce the linear operators C and C 1 describing a perturbation of the stress field inside 
and outside the inclusion Vi (x £ R d ) 



The right-hand side of Eqs. (3.17) and (3.18) can be also estimated in the spirit of subtraction technique 
(3.2) according to the next scheme 



for x e Vi and x ^ Vi, respectively; here x TO = argmin y |x - y (y G i'„ x ^ v{), C = C° , g = <T;, T7, 
and the tensor Q^x) is defined analogously to Eq. (3.3). 

We constructed the solution (3.16) for a perturbation of the stress field inside and outside the inclusion 
Vi in the operator form obtained by the method of volume integral equation (VIE) for an arbitrary effective 
field <Xi(x), (x G Vi). However, this operator could be created by any another numerical method such as, 
e.g. the finite element analysis (FEA). The main difficulty in such a case is a generation of prescribed 
effective field <Xi(x), (x 6 «,) which will be considered in the next subsection. 

3.2 Creation of prescribed stresses by the FEA 

Construction of the operator C a (3.16) anticipates a creation of prescribed effective field Wi(x), (x € 
Vi) in the absence of the inclusion Vi. In the case of the FEA employment, it can be done by prescribing of 
either some boundary condition at the boundary of a large sample or some eigenstress inside this sample. 
We will only consider the second way. 

The problem is to find a fictitious /3i(x) generating a prescribed stress <r(x) in an arbitrary fictitious 
ellipsoidal inclusion x € vq (which has no connection with the correlation hole «?■) with the elastic 
modulus M(x) = M*- ' and an indicator function Vq. 




(3.16) 



(3.17) 



(3.18) 




(3.19) 



(3.20) 




(3.21) 



The equation (|3.2ip can be recast in the form 




(3.22) 



On the thermo-elastostatics of heterogeneous materials. II 11 

Except for notations, the Fredholm integral equation of the second kind (|3.22[) coincides with the direct 
equation for estimation of stresses <x(x) produced by the field <f(x) inside the heterogeneity (x G vq) 

<t(x) - Bo^(x) + B J r(x - y)[M 1 (y) < T(y) - M 1 (x) < r(x)]Vb(y)dy (3.23) 

Indeed, Eq. (|3.22|) is reduced to Eq. (|3.23|) (with /3 ± = 0) by the replacement of the notations: (3 1 — » Mi<r, 
-Qo 1 ^ -> MiBoo 7 , Q lr -> MiB r, where B = B (v ) is defined analogously to (2.14). 

In the case of the successive approximations method, we need to evaluate the right-hand side of the 
equation 

M n+11 (x) = -Qo V(x) + q 1 J r(x - y )[M nl (y) - M nl (x)]v (y)dy (3-24) 

with the usual use of the driving term as an initial approximation 

M 0] (x) = -Qo^(x). (3.25) 

The volume integral equation (|3.2ip is reduced to the regular representation, which has no singularities 
and can be also presented in the form adopted for using of the FEA 

/3[" +1 l(x) = Q - V(x) + M nl (x) - Qo-^/^x) (3.26) 

where the operator 

-R*M nl (x) = /r(x - y)M nl (y)^o(y)dy (3.27) 

presents the stresses produced by the intermediate eigenstresses /^"'(x) in the inclusion x G vq (see 
Eq. (13.2111 ). Obviously, the mentioned stresses can be easy estimated by the FEA. Thus, an operator 
representation of the solution (|3.2ip 

/3(x) = r-Wtx) ee J r- x (x - y)W(y)V (y)dy (3.28) 

can be considered as found by the FEA (or by any other numerical method providing a solution of the 
regular integral Eq. (|3.22p . The next step for FEA utilization is obvious. We introduce the real inclusion 
Vi into the fictitious ellipsoid vo (such that all desirable area for the stress estimation is placed inside vo) 
and estimate the real stresses <r(x) which can be considered as found in Eq. (3.16). 

The VIE and FEA methods have a series of advantages and disadvantages (considered, e.g., in Rcf. 
[6]), and it is crucial for the analyst to be aware of their range of applications. 

4. Estimation of both the effective field and effective elastic moduli 

The new general integral equation (2.1) can be rewritten in terms of the operator representation C a 
(3.16) 

<r(x) = <<x)(x) + J [.CfaXx) - (Cfr))M]dy (4.1) 

while conditional averaging of Eqs. (2.72) and (2.73) leads to the following representation for the mean 
of the effective field in the fixed inhomogeneity x G Ui 

(5%(x) = <<x)(x)+ f [£J((5l;^,x i > ? )(x)^„x,|; I , i ,x i ) - £J((<f),)(x)n<«>(x,)]dx,, (4.2) 
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where (<f|;«,,Xj} 9 = (<x|; v i} x i ) (? (y) is a conditional statistical average of <r(y) varying along the fixed 
heterogeneity y <E v q at the fixed i>i while (<x) 9 = (<r) g (y) is a statistical average <f(y) inside the hetero- 
geneity y 6 v q . No confusion will arise hereafter in definition of the operator T> (CZ, C q ) with the kernel 
2>(x, y) on the inhomogeneous functions g(x) (e.g., g(y) = (W) q (y), y £ V k ;V k = V h V q ,Vi + V q ) 

2>(g)(x) = /l>(x,y)g(y)\4(y)dy. (4.3) 

The operator X> is reduced to the tensor D(x) on the constant functions g(x) = g =const (x £ Vk) 

2>(g)(x) = D(x)g, D(x)= f -D(x,y)V k (y)dy. (4.4) 

The integral in the right-hand side of Eq. (4.2) can be decomposed as 

(5%(x) = (er)(x) + J, J = Ji + J 2 + J 3 , where (4.5) 
Ji = J [^g((5 r |;«i,x i ) g )(x) -£^((CT) g )(x)]^(u 9 ,Xg|;u i ,x i )dx g , (4.6) 

J 2 = | C^(W) q )^)[^v qiXq \;v llXl ) -nW(x,)][l - ^ g (x g )]dx 9 , (4.7) 

J 3 = - / £J((ff),)(x)nW(x,)^(x,)dx,. (4.8) 

The absolutely convergent integral in Eq. (4.2) is decomposed in Eq. (4.5) just for subsequent presentation 
obviousness; because of this, the absolute convergences of integrals Ji, J 2 and J3 are not considered. 

In the framework of the quasi-crystalline approximation (2.22) (y G v q ) 

(<rU,x 4 ) 9 (y) = (<r) q (y), (4.9) 
<if|;« ijXi >,(y) = (W) q (y), (4.10) 

and Eq. (4.5) is simplified (J = J2 + J3) 

(^(x) = (<t)(x) + I C a q ((W) q )( X )y(v q ,K q \;v t ,^) - n^( Xq )]d Xq , (4.11) 

For statistically homogeneous media when n(x q ) = wS q > =const., it is logical to assume the acceptance of 
the additional hypothesis of "ellipsoidal symmetry" (2.24), which leads to the same simplification as in 
Subsection 2.4 (J 2 = 0) 



C^(W} q )( X )y(v q , Xq \; Vl ,^) - nM(x,)][l - t£(x,)]dx, = 0, (4.12) 
which can be presented in an equivalent form exploiting Green's function 

y - r(x-y)[(M lCT ) 9 (y)+/3 1 (y)]bK,x ? |; U , ) xO-n^]y g (y)[l-^°(y)]dy = 0. (4.13) 
Then Eq. (4.5) leads to (J = J 3 ) 

(5%(x) = <<t)(x) - / n(«)£J((ff) 9 )(x)^ (x,)dx,. (4.14) 
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The accuracies of the assumptions (4.9), (4.10) and (4.12) will be estimated in Section 6. However, we will 
perform subsequent solution of Eq. (4.11) rather than Eq. (4.14) by keeping in mind that the hypothesis 
H3 (2.24) is accepted. 

Obviously, the regular integral equation (4.11) has no singularities and can be solved by a direct 
quadrature method with formal representation of the solution 

(ff) i (x)=T 1 -*< ff )(x) (4.15) 

Although the direct quadrature method usually causes no problems of accuracy, for a large number 
of unknown variables N its 0(N 3 ) cost dependence can lead to surprisingly long computing time. The 
obvious way of reducing this cost is to construct an iterative scheme which will be considered now. Namely, 
Eq. (4.11) will be solved by the iteration method when the initial constant effective stress (<r) (x) = 
(<r) 9 =const. is estimated from the classical approach (2.18) and (2.24). Indeed, for (cr) =const., the 
statistical average of stresses inside the inhomogeneity v q is found to be inhomogeneous [see Eq. (2.12)] 

(«7> 9 (x) = B,(x)(ff), + C,(x) ) (r,> ? (x)=R*(x)(<f) 3 + F*(x), (4.16) 

in a general case of inhomogeneity of v q ; here = v q 1 R q , F q = v^Fg. In such a case, the right-hand 
side of Eq. (4.11) generates inhomogeneous field (er) i (x). For elimination of this difficulty, we will use the 
additional condition of the effective field hypothesis Hlb (2.10), when a perturbation introduced by the 
inhomogeneity v q is defined by the strain polarization tensor (r/) averaged over the volume v q (x £ R d ) 

£*«cr> fl )(x) = v q T q ( X - x q )(r,) q , (4.17) 

Then Eq. (4.11) in the framework of the hypothesis H3 (2.24) is reduced to the classical representation 
for the effective field (<r)-(x) 

(ff>.(x) = (ir> + Q 1 (x)X; C (9)^) 9 , (4.18) 

9=1 

which is homogeneous just for an additional assumption of an ellipsoidal shape of the excluded volume 
u?. Combining the averaged Eqs. (4.16) and (4.18) leads to the final representations for the averaged 
tensors of both the effective field and stain polarization 

<7x)[ 01 = ( < T) + Ql ) [I-(R u QV)]" 1 ((R u y)<<T) + (F 1 V)), (4.19) 
(x) - B,(x)(«t) + C 1 (x)+B,(x)Q»[I-(R"QV)]- 1 ((RT)(<t) + (FT)), (4.20) 
(77>f ] (x) = R,(x)(«t)+F,(x)+R ! (x)Q»[I-(R"QV)]- 1 ((RT)(<t) + (FT)), (4.21) 

which will be considered as the initial approximation of the next equations (x S Wj, y£« ? ) 

(<f)|" +1] (x) = (*) + J ^«^ n] )(x)b(t; fl ,x,|;« i ,x i )-nW(x,)]dx„ (4.22) 
<r?>i n+1I (y) =7V^}' n+11 (y)+F^y), (4.23) 

where lZ q = Mi (I + C q ) and Eq. (4.22) in the case of the assumption (4.12) is reduced to the following 
one 

<7x)f +11 (x) = <">-/ ^'((»7)J ll )(x)nW(x 9 )^(x 9 )dx g . (4.24) 
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The system (4.22) and (4.23) can be formally presented in an operator form (tj)(x) = (<t) + /C((t/))(x) 
(the indexes are dropped for simplicity). It suggests the Neumann series form for the solution 77 of (4.22) 
and (4.23) [compare with the solution (4.15)] 

fa>,(x) = lim faM^x) = R*(x)<<r) + F*(x), (4.25) 

n — >oc 

which yields the final representations for the effective properties 

M* = M (0) + (R*V), /3* = /3 (0) + (F*V). (4.26) 

A convergence of the sequence (^["^(x) (4.25) is analyzed analogously to the sequence (3.14). 

In Eq. (4.26) we used an obvious connection between the phase average (gV) (g = (7,6,77) and the 
averages inside the representative inclusions Vk £ (h = 1, . . . , N) 

N 

<g^H5>( fc >(g> fc , (4.27) 
fc=i 

which is only fulfilled for statistically homogeneous media subjected to the homogeneous boundary con- 
ditions. If any of these conditions is broken then it is necessary to consider a generalization of Eq. (4.27) 
in the form of Eq. (3.291). However, the mentioned class of nonlocal problems is beyond the scope of the 
present paper. 

5. Qualitative analysis of some basic hypotheses and propositions 

5.1 Analysis of the proposition 1 and hypothesis HI 

The hypothesis HI is widely used (explicitly or implicitly) for the majority of the methods of 
micromechanics even if the term "effective field hypothesis" is not indicated. For example, Buryachcnko [6] 
demonstrated that hypothesis HI is exploited in the effective medium method, generalized self-consistent 
method, differential methods, Mori-Tanaka method, the MEFM, conditional moments method, variational 
methods, and others. These are a lot of other methods using the hypothesis HI differ one from one another 
by some additional specific assumptions. 

It should be mentioned, that the domain of the operator £^((t7)|^)(x) (3.18) is a whole space x £ R d , 
and, because of this, some points of the area x £ Vi in Eq. (4.22) can be uncovered by the heterogeneities v q 
and, therefore, the effective stress (<x) (x) (4.22) will depend on the stress perturbations C q ((T/)q™')(x) 
in the vicinity x £ Vi \ v q of the area v q rather than only on stress distributions in the inhomogeneity 
v q and max x |x — x q | = 3a (x £ Vi \ v q ) for the identical spherical inhomogencities of the radius a with 
an isotropic statistically distribution of their centers. Thus, we obtain a fundamental conclusion that 
effective moduli in general depend not only on the stress distribution inside the inhomogencities but also 
on the stresses in the vicinities of inhomogencities (compare with the proposition 1). However, if our 
estimations utilize Eq. (2.2) containing only average strain polarization r\ q [rather than (r(x — y)r](y)) q 
in Eq. (4.1)] as a renormalizing item then an influence of stresses in the vicinities of inhomogencities 
is degenerated. At the same time, using Eq. (4.1) leads to the necessity of evaluation of stresses in the 
inhomogeneity vicinities even for a statistically homogeneous field of ellipsoidal homogeneous inclusions 
(it will be quantitatively demonstrated in Section 6). Moreover, a fundamental deficiency of Eq. (2.2) is 
the dependence of the renormalizing item T(x — x q ){rj} only on the average stress polarization tensor 
(77) while a corresponding item (T(x — y)rj(y)) q (x q ) (y £ x 9 ) in the new Eq. (4.1) explicitly depends 
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on details distribution (r]\v q ,x q )(y) (y G x q ). Because of this, the averaging methods used in Eq. (2.2) 
as a starting element conserve the mentioned deficiency of Eq. (2.2) (at least in some elements of these 
methods). For example, Chen and Acrivos [36] have estimated the effective elastic moduli through the 
accurate evaluation of binary interactions of inclusions without the hypothesis HI, c.i. in our notations 
the operator C q ((tf|; Vi, Xj) g )(x) (4.2) was estimated at the condition (<7j, g ) = (er) [compare with Eq. 
(2.17)]. However, the operator ((<f))(x) (4.2) was estimated at the approximation (4.17) that implicitly 
implies the use of both the hypothesis HI and Eq. (2.2) (see for details Subsection 10.2.2 in [6]). 

On the other hand, although the method (4.18)-(4.21) allows the inhomogeneous statistically aver- 
aged tensor {rj) q (y) (y S v q ) (4.21), but Eq. (4.21) containing the item (R^Q ^) generated by the second 
summand in the left-hand side of Eq. (4.18) depending only on the average stress polarization tensor (77). 
As a consequence of this the final classical representations of the effective properties (2.19), (2.20) and 
(2.25) depend only on average stress concentrator factors Ri and while the effective properties (4.26) 
explicitly depend on the inhomogeneous tensors Ri(x) and F^(x) as well as on detailed distribution 

(v) q (y) (y g v q ) (4.23). 

Moreover, the detected explicit dependence of the effective properties (4.26) on the detailed stress 
concentrator factors R;(x) and F^(x) rather than on the average values R^ and F^ allows for an aban- 
donment of the hypothesis Hlb [or (4.17)] whose accuracy is questionable for the inhomogeneous (e.g., 
coated) inclusions. In such a case the statistical average effective field estimated by Eq. (4.22) is found 
to be inhomogeneous that discards the hypothesis Hla. Quantitative estimations of the result of this 
abandonment of the hypothesis HI will be performed in Section 6 in the framework of the hypothesis 
H3 for some particular cases of fiber composites. 

5.2 Analysis of the proposition 2 and hypotheses Hlb and H3 

As was mentioned, forfeiting of the effective filed hypothesis Hla by the additional perturbation 
hypothesis (4.17) leads to Eq. (4.18) providing homogeneity of the effective field estimation for the 
ellipsoidal excluded volume Moreover this estimation of the effective field (and, therefore, of effective 
moduli) is invariant with respect to the size of the ellipsoidal excluded volume v®. However, the additional 
hypothesis (4.17) is exactly fulfilled only for the homogeneous ellipsoidal inhomogencity v^. For both the 
inhomogeneous and noncllipsoidal inclusions the equality (4.17) is just an approximation and the new 
general equation (4.1) has an advantage with respect to the popular one (2.2) only based on average 
strain polarization r\ q [rather than (r(x — y)r;(y)) ]. Then the size of the excluded volume i>° will impact 
on the effective field (4.22). Indeed, if the radius of the excluded volume in Fig. 1 increases from 2a 
to 3a then the long distance of influence zone of the inhomogeneity v q on the effective field (ct)-(x) will 
increase from the value x — x g | = 3a (as in Fig. 1) till |x — x q \ = 4a. This influence will be quantitatively 
estimated in the next section. 

A popular explanation of acceptance of the "ellipsoidal symmetry" hypothesis (2.24) is that this 
hypothesis just simplifies Eq. (2.23) reducing this equating to Eq. (2.25) which does not contain the 
integrals. In a similar manner, a destination of the assumption of the ellipsoidal shape of the excluded 
volume vf in the hypothesis H3 is that this hypothesis just simplifies Eq. (4.18) by the use of analytical 
known tensor Q° (expressed through the Eshelby tensor (3.4)) which is exploited instead of a general 
tensor Q°(x) found numerically (see e.g. Subsection 4.7.4 in the book [6]). However, the both mentioned 
assumptions of the hypothesis H3 have a fundamental conceptual sense rather than only an analytical 
solution of some particular problem. Exploiting the Eshelby tensor concept in Eq. (4.18) (and in the 
MEFM) is based on the ellipsoidal shape of the correlation hole w° rather than on the inclusion shape 
Vi. An abandonment of either the assumption of the vf's ellipsoidal shape or "ellipsoidal symmetry" 
hypothesis (2.24) with necessarily leads to the inhomogencity of the effective field Wi acting on the 
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inclusion Vi that is prohibited for the classical version of the MEFM. However, Buryachenko [6] (see 
Section 9.4) proposed a method for solution of Eq. (4.18) based on the general integral equation (2.2). 
Namely, acting on Eq. (4.18) by the operator & = (I TMi) -1 ^ yields 




Fig. 1: Schematic mutual 
placements of v q , Vi and v® 

N 

(fj),(x) = Ri(x)(*) + F 4 (x) + R Q0 (x) 5>,) C <*>. (5.1) 

9=1 

where the numerical estimation scheme of the tensor R^°(x) for any shape of v® was considered in 
Subsection 9.4.4 in Ref. [6]. Volume averaging of Eq. (5.1) over the heterogeneity Vi and summation 
over the inclusion number i lead to the average strain polarization tensor (rj) and, therefore, gives the 
representations for the effective properties 

M* = M (0) + [I- (R Q0 nOr l (R 1 T) ) (5.2) 
f3* = /3 (0) + [l-(R QOv V)]~ 1 (F v V). (5.3) 

If (and only if) the correlation hole vf is chosen as an ellipsoid (homothctical, for example, to Vi) then 
the tensor Q?(x) = Q° = const., R 4 Q0 (x) = v l (w ? °)- 1 R. i (x)Q^, B 4 Q0 (x) = B°(x)Q°, and, therefore, Eqs. 

(5.2) and (5.3) are reduced to the known representations (2.20) and (2.25) with the constant tensor Q° 
depending on the orientation of the correlation hole v®. An advantage of Eqs. (4.22)-(4.26) with respect 
to Eqs. (5.2) and (5.3) are that Eqs. (5.2) and (5.3) are fundamentally limited by analysis of statistically 
homogeneous media subjected to the homogeneous boundary conditions while the system (4.22) and 
(4.23) can be easily generalized to the statistically inhomogeneous media. Moreover, the method (5.1)- 

(5.3) conserves a fundamental deficiency of the general integral Eq. (2.2) containing the renormalizing 
item defining only by the averaged strain polarization tensor (rj). Estimation of the tensor R^°(x) 
in an auxiliary model problem with homogeneous fictitious eigenstrain in v® implies that influence of 
surrounding inclusions v q (x g S vf) is defined by the average strain polarization tensors (rj) rather than 
its detailed distribution (»/) (x) (x G v q ). Impact of the last improvement on the estimated effective 
properties will be considered in Section 6 for the circle vf although, of course, analysis of any shape of vf 
present no additional difficulties for the new method (4.22)-(4.23) as opposed to the method (5.1)-(5.3) 
requiring evaluation of the supplementary tensor R^°(x). 

5.3 Analyses of the hypotheses H2a and H2b 

As it was noted in Subsection 2.3, the hypotheses H2a and H2b are not conceptually dependent 
on the hypothesis HI and can be applied in general case even if the hypothesis HI is violated. Indeed, 
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at obtaining of Eq. (2.18) we already used the hypothesis HI because Eq. (2.18) contains the objects 
Rg(tf) 9 + F q instead of their operator generalization jC^(rj)(x) (3.18) which does not use the hypothesis 
HI. However, even in this case the effective field (W) i (x) (2.18) is an inhomogeneous function of the 
coordinate x e Uj. In actual truss, a subsequent averaging of Eq. (2.18) over the inclusion Vi is tantamount 
to a secondary using of the hypothesis HI that is not necessary and can be avoided. However, such an 
inhomogeneity of (er^x) (2.18) is beyond the scope of the current study and will be analyzed in other 
publications. It is correctly noted that application to Eq. (2.18) of the simplified hypothesis H2b [with 
= ISij (2.22)] in an accompany with the hypothesis H3 docs not fall outside the scope the hypothesis 
HI and reduces Eq. (2.18) to Eq. (4.18) with subsequent obtaining of the known representations for the 
effective properties (2.20) and (2.25). However, the eventual abandonment of the hypothesis HI can be 
done before the use of the hypotheses either H2a or H2b as it was performed in Eq. (4.2). Then the 
following solution of Eq. (4.2) by the use of the hypotheses either H2a, H2b or H3 does not lead to 
the necessity of using the hypothesis HI that will be quantatively demonstrated in Section 6 at some 
numerical examples. 

6. Numerical results 

With the non-essential restriction on space dimensionality d and the shape of inhomogencities we will 
consider 2-D problems for composites reinforced by cylindrical infinite fibers. The domains of inclusions 
Vi are discretized along the polar angle and the radius in the local polar coordinate system with the 
centers x*. Then the points 

\(r,<p) | (p-1 — <<p<p-r- , {q-1 - <r<q-\ (6.1) 
t I l mm) 

[p = 1, 2 . . . , I; q = 1, 2, . . . , m) represent the elements of T? q of the meshes fii (i = 2, . . . , n) that is not 
optimized, but is efficient. Moreover, the square meshes 

{(*!,x 2 ) T | (p-l)^. <Xl<p ^, {q -l)^ <X2<q ^y (6.2) 

where x\, xi are local coordinates with origins at the fiber centers, will be used for stress estimation inside 
and outside the fiber. We will use piecewise-constant elements of the meshes which are not very cost- 
efficient but are very easy for computer programming, and the discretization (6.2) permits the analysis 
of nonregular inclusion shapes. For simplicity estimation of integrals involved we will utilize the basic 
numerical integrations formulas of Simpson's rule and trapezoidal rule for the uniform (6.1)-(6.2) and 
nonuniform meshes considered below, respectively. 

We detected that in the concrete examples of high matrix-inclusion elastic contrast considered and 
some others, the standard popular iterative schemes (3.14) may diverge or converge very slowly (i.e. the 
iteration scheme (3.12) does not work in general) so that an implementation of the improved algorithm 
proposed in this paper becomes more complicated. In such a case, following Refs. [37], [38], we introduced 
the subsidiary grid of the support points Cj at the centers of each elements additionally to the nodal 
points Sj at the apexes of elements. After determining in this way at all the points £j the values of the 
function rjN(£j), we find its values at the nodal points rf^(sj) by linear interpolation, and so on (see 
details in Rcf. [38]). Moreover, instead of the point Jacobi iteration method displayed in Eq. (3.14) we 
use the accelerated Liebmann method (called also extrapolated Gauss-Seidel method) which is usually 
"faster" than the point Jacobi method, and has the computational advantage that it does not require 
the simultaneous storage of the two iterations Vtk+i) an d V(k) ( see i e -S-i [39]). The convergence of the 
scheme (3.14) is provided by their modification r]^- k+1 ^ = ^[r)^ + + /CsT/W] (see for details Ref. [38]). 
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It should be mentioned that in forthcoming numerical examples we will use only the iteration scheme 
described above. Comparative analysis of this scheme with other known iteration schemes is beyond the 
scope of the current paper. Moreover, although the convergence of this method was rigorously proved in 
Ref. [38] for the elastic problems of an arbitrary dimensions, we will demonstrate its effectiveness only 
for 2-D problems; the analysis of 3-D problems is beyond the scope of this paper. 

We consider a pure mechanical problem (/3 = 0) and assume the matrix is epoxy resin (l/ -* = 
(3fc(°\ 2/LtW), fcW = 3.83 GPa and ^ = 1.27 GPa) which contains identical circular glass fibers (L (1) = 
(3fcW,2^W), fcW = 34.3 GPa and ^ = 31.3 GPa). If the pair distribution function 5 (x ?; - x TO ) = 
<p(vi, Xj|; v m , x m )/inS k > depends on x m — x;| it is called the radial distribution function (RDF). Two 
alternative RDFs of inclusion will be examined (see Refs. [40], [41]) 

g(xi - x 9 ) = cp(vi,Xi\;v q ,x q )/n^ = H(r - 2a), (6.3) 
g(xi - xj) = H(r - 2a)|l + ^ [tt - 2siir 1 (^) - I-yJl-J^H(4a - r)} (6.4) 

where H denotes the Heaviside step function, r = |xj — x 9 is the distance between the nonintersecting 
inclusions Vi and v q , and c is the volume fraction of fibers of the radius a. The formula (6.4) takes into 
account a neighboring order in the distribution of the inclusions. 

At first wc will perform our evaluations for composites with homogeneous fibers described by the RDF 
(6.3) in the framework of the hypotheses H2b and H3. Influence of the effective field hypothesis HI is 
considered by comparison of statistical averages of stresses in the fibers estimated by the classical approach 
(cr)° ld (x) =const. (4.15), (4.18) as well as by the proposed one (<x)" ow (x) (4.22), (4.23). Wc considered 

a volume fraction of fibers c = 0.65 and evaluated the stress perturbations C q ((rj} q n ^)(x.) (4.22) in the 
vicinity {x| max x |x — x g | = 3a} of the area v q rather than only a stress distributions in the inhomogeneity 
x G v q . Then (<x)° ld and (cr)" ew (x) differ from one another no more than 0.09% that coincides with a 
computational error realized in the method (4.22), (4.23) for two different meshes (6.2) with I = 15 
and I = 30. Thus, we qualitatively proved that in the considered example both methods the old (4.15), 
(4.18) and new (4.22), (4.23) ones which are based on the classical (2.2) and new (4.1) general integral 
equations, respectively, lead to the same numerical results. This conclusion quantitatively confirms the 
Proposition 1) establishing an equivalentness of Eqs. (2.2) and (4.1) for statistically homogeneous fields 
of homogeneous ellipsoidal heterogeneities subjected to the homogeneous boundary conditions. Now we 
will consider an influence of incorrect using of the operator C q ((rj) q n] )(x) (4.22) when only x G v q are 
considered, which means that stress perturbations introduced by the moving inhomogeneity v q in their 
vicinity {x| a < |x — x. q \ < 3a} are neglected. For this purpose, the means of stress concentrator factors 
(0 = 0) 

( < T) i (x) = B*(x)( CT ), (6.5) 

defined analogously to Eq. (4.25), will be estimated. In Fig. 2 the components -B^iifx) demonstrating 
maximum dependence on x = (xi,0) T arc presented for the initial [i?22ii( x ) — ^2211 =const.], second 

[-^2211 ( X )L forth [^2211 ( X )L an d tenth [-B^n^x)] iterations of stress concentrator factor. A fast conver- 
gence of the proposed iteration method can be seen: the tenth iteration differs from the ninth, fourth, 
and initial approximations by 0.013%, 0.084%, and 2.9%, respectively. In so doing, the difference 2.9% 
essentially exceeds the possible errors of both the calculations and iteration scheme. Thus, for statistically 
homogeneous fields of homogeneous circle inclusions subjected to the homogeneous boundary conditions, 
the old and new approaches based on the backgrounds in the form of Eqs. (2.2) and (4.1), respectively, 
lead to equivalent results. Thus, in the case of the background (4.1) we must estimate the stress pertur- 
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bation in the vicinity {x| a < 



to the proposition 1 obtained at the use of the old background (2.2). 



< 3a} of the moving inhomogeneity v q . This statement contradicts 
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Fig. 2:^2211 (xi) vs Xi/a in homogeneous 
fiber: curves 1, 2, 3, 4 for fc = 0, 2, 4, 12, 
respectively 

We are expected to get a larger difference of the backgrounds (2.2) and (4.1) for composites reinforced 
by either noncllipsoidal or inhomogencous inclusions demonstrating essentially inhomogencous stress dis- 
tribution inside inclusions even in the framework of the hypothesis HI. The interphase is usually the 
product of processing conditions involved in composite manufacture. In relation to this problem, an appli- 
cation of the concept of functionally graded materials by Hirai et al. [42] for description of the interphase 
whose moduli may vary continuously is worthy of notice. Along this line one may, for instance, refer to 
the works [43-47] concerned with the spatially nonuniform properties of interphase. Just for concreteness, 
we assume that fibers contain the cores of the radius a c < a with the constant moduli L/ 1 - 1 =const while 



the moduli L m (x) in the interphase with the coating thickness h 



the radial distance 



L int (r) =L (o) + (L (D _ L (°))(a-r)// l 



a — a c are taken to vary linearly with 

(6.6) 



For demonstration of maximum difference between the old and new approaches, we will consider in detail 
a thick coating with the relative coating thickness h/a c = 0.5 although other ratios h/a c will be also 
analyzed in a few comparative examples. At first, we will analyze results obtained in the framework the 
hypotheses H2b and H3 for the RDF (6.3). In Fig. 3 the iterations S^fx) (k = 0, 2, 4, 12) at the axis 
x = (xi,0) T are presented for c = Tra 2 n = 0.65. The initial approximation ^[^(x) corresponding to 
the classical estimation (4.20) and using the old background (2.2) reveals their essential inhomogeneity 
(14%) even in the framework of the effective field hypothesis HI. The new background (4.1) allow the 
use of this inhomogeneity for refinement of the rcnormalizing item in Eq. (4.14) without exploiting of 
the hypothesis HI. The twelfth iteration S*|]j(x) differs from the initial approximation ^[^(x) by 
10.7% while the 12th and 11th iterations are distinguished from one another by 0.9%. Of even greater 
difference of results obtained for the backgrounds (2.2) and (4.1) is observed for the component B*^ 2 (' x ) 



(k = 0, 2, 4, 12) at x = (x x , 0) T in Fig. 4. Indeed, S*|a J 2 (x) > at any x = (x 1} 0) T while B*!$ 2 (x.) < at 
—0.75a < x\ < 0.75a. Again, the proposed iteration method converges rapidly and wc contend that 12th 
iteration provides a difference from 11th iteration of 0.2%. 
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Fig. 3: B^[~^ 1 (xi) vs X\ja in inhomogeneous Fig. 4: B^2 2 {xi) vs x l/ a in inhomogeneous 

fiber: curves 1, 2, 3, 4 for fc = 0, 2, 4, 12, fiber: curves 1, 2, 3, 4 for fc = 0, 2, 4, 12, 

respectiveiy. respectiveiy. 

So much prominent and systematic differences of the old and new approaches are based on the 
abandonment from effective field hypothesis HI in the new approach. We estimated a tensor of effective 
stress concentrator factor 

<<x}(x)=B*(x)<<x}, (x=( a;i ,0) T ) (6.7) 
and presented the components of their fc-th approximations i? lln (x) and -B 1122 (x) (fc = 0,2,4,12) in 

*[12] *[0] *[H] 

Figs. 5 and 6, respectively. £? lm (x) differs from both the classical £? lm (x) (4.15), (6.7) and i? ml (x) 

jjj h2] 

on 8.1% and 0.03%, respectively, while -B lm (x) varies along x\ over 2.2%. However, we can observe in 
Fig. 6 a significantly more dramatic situation with the component £?ii22( x ) where all iterations differ by 
a sign from the classical one -B,i22( x ) (4.15), (6.7) almost at all values \x\\ < a. 

We now turn our attention to the analysis of the size of the circle excluded volume v® with the radius 
a° on the stress concentrator factor B*^ fe '(x) also for the radial distribution function (6.3) reducing Eq. 
(4.22) to Eq. (4.24). We will compare the estimation of B*^(x) carried out for a = 3a with previously 
obtained results for a = 2a (see Figs. 3 and 4). The components .B^^'-^x) and -Eq^M ( x = (• T ij0) T ) 
are presented in Figs. 7 and 8, respectively, for fc = (curves 1) and fc = 12 for both a = 2a (curves 
2) and a° = 3a (curves 3). The RDF (6.3) provides the "ellipsoidal symmetry" hypothesis H3 (2.24) 
and, because of this, the classical representations for B*'°'(x) is invariant to the size of w° while S^[\^(x) 
and Sij^fW estimated by the new approach (4.24) for a = 2a and a = 3a differ at x\ = one from 
another by 3.1% and 50%, respectively. Finally, we compare the influence of the RDF (6.3) and (6.4) at 
a = 2a on estimation of B*' fc '(x). Needless to mention that B*'°'(x) (4.20) is invariant with respect to 
the RDF while B*} 1 ^ 1 (x) and B$${x) estimated for the RDF (6.3) (curve 2) and (6.4) (curves 4) are 
distinguished by 3.7% and 33%, respectively. The indicated differences demonstrating fundamentally new 
effects inherent in the new approach (4.22)-(4.25) far exceed the iteration error between 11th and 12th 
iterations which are less than 0.03%. 
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Fig. 5: Effective field concentrator factors 

-B 1111 (.Ti) vs xi/a in inhomogeneous fiber: 
curves 1, 2, 3, 4 for fc = 0, 2, 4, 12, respectively. 
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B 112 2( x i) vs X\/a in inhomogeneous fiber: 
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Fig. 8: Fig. 8. 5*122 (^l) vs x i/ a for the 
different a and RDF: curves 1 (k = 0), 2 
- 2a, fc = 12], 3 [RDF (6.3), 



1.0 



[RDF (6.3), a 



Fig. 7: Bnii(xi) vs xi/a for the different a 
and RDF: curves 1 (k = 0), 2 [RDF (6.3), 
a = 2a, & = 12], 3 [RDF (6.3), 

a = 3a, fc = 12], 4 [RDF (6.4), a = 3a, k 

a = 2a, k = 12]. a = 

Just for completeness, we will estimate an influence of the interphase thickness h (6.6) on the stress 
concentrator factor B* [0] (x) and B* [12] (x) (x = (xi,0) T ) for the RDF (6.3) with a° = 2a. In addition to 
Figs. 3 and 4 displaying the results for h/a c = 0.5, we are demonstrating the similar estimations B*^ 2 (' x ) 



= 12], 4 [RDF (6.4), 
2a, k = 121 



(k = 0, 12) for h/a c = 0.1, 0.25, and h/a c = 1 in Fig. 9. As can be seen, B^ 2 \{a) and B*^ 22 (a) differ one 
from another by 20.5%, 30.3%, and 33.1% for h/a c = 0.1, 0.25, and h/a c = 1, respectively. The similar 



,*[12] 
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differences for the components B*[°\(a) and (a) are 6.0%, 7.2%, and 8.5%, respectively. 

We complete our numerical analysis by estimation of isotropic effective moduli L* = 2fc* 2 |Ni+2/i* 2 jN2 



"[2]) 
,(0) 



(Ni = S (g) 5/2, N2 = I Ni). For the fiber composites it is the plane-strain bulk modulus k^ (and 

j) - instead of the 3-D bulk modulus k$ - that plays the significant role: feS = fcS + /%] /3j 

/.tj 2 j' = Mpi • A i *// i< "°' ) are presented in Fig. 10 for both the classical approach [corresponding to the stress 

concentrator factors B*'°'(x)] and new one [corresponding to the 12th iteration B*^ 12 '(x)]. As can be 
seen, the distinctions between two approaches equal 3.8% and 12.0% for c = 0.65 for fx*/fj,^ and k*/k^°\ 
respectively. In so doing, the stress concentrator factors in these approaches at the point x\ = a of fibers 
can differ on 30% and, moreover, these estimations for the different approaches can have the different 
signs at other same domains of vi (see Fig. 4). Thus, stress concentrator factors are significantly more 
sensitive values to the choice of the approach than effective elastic moduli. 
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7. Conclusion 



We have proposed the new background of micromechanics based on the new general integral equation 
(4.1) which does not use the central concept of classical micromechanics such as effective field hypothesis 
HI. The eventual abandonment from hypothesis HI has made a rejection of the satellite hypothesis 
H3 possible. If statistical averages of stresses in the heterogeneities can be considered as homogeneous 
ones then the new approach is degenerated into the classical approach (4.19)-(4.21). However such an 
assumption is approximately appropriate only for statistically homogeneous fields of homogeneous ellip- 
soidal inhomogencitics subjected to homogeneous boundary conditions and fulfilled at the conditions of 
quasi-crystallitc approximation 3. If any of the indicated conditions is broken then an appearing inho- 
mogeneity of stress fields in the inclusions lead to one of two possible sources of inhomogencitics of the 
effective field which, in turn, generates an additional inhomogencity of stress fields inside inclusions and 
so on. For example, if all above-listed conditions are satisfied but the closing hypothesis H2b is replaced 
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by the hypothesis H2a taking the binary interaction of inclusions into account then the first sort of in- 
homogeneity of the effective field (<f ) (x) is generated by the binary interaction of inclusions [see the item 
C? q ((<x |; Vi, Xi)^)(x) in Eq. (4.6)] even if this interaction is approximately estimated through the matrix 
Z as in Eq. (2.18). However, even in the framework of hypothesis H2b, the replacement of ellipsoidal 
homogeneous inclusions by either the nonellipsoidal homogeneous ones or inhomogeneous (e.g. coated) 
ellipsoidal inclusions with necessity leads to the second sort of effective field inhomogeneity produced 
by the fundamentally new renormalizing item (C* (r))) (x.) (4.1). This new renormalizing item is directly 
dependent [in opposite to the classical Eq. (2.2)] on inhomogeneity of stress fields inside the inclusions 
that has lead to detection of fundamentally new effects in micromechanics such as dependence of stress 
concentrator factors estimated (see Figs. 7 and 8) on both the RDF and size of the excluded volume even 
in the framework of hypotheses H2b and H3. 

The modeling and simulation of random nano- and microstructures are becoming more and more 
ambitious due to the advances in modern computer software and hardware that is stimulated by a real 
challenge of modern material science and technology. The researches can forget about restrictions of 
analytical solutions (such as, e.g., Eshelby tensor and hypothesis HI) and use the numerical solutions 
which they need. It is expected to get all the more differences between the old and new approaches 
than inhomogeneity of the stress concentrator factors B^(x) (x £ m) would be larger. So, for the square 
inclusion with the smoothed vertexes and the finite cylindrical fiber the components Bj(x) can vary by 
factors of four and ten (compare with Fig. 2), respectively (see Subsections 4.2.4 and 18.3.2, respectively, 
in Rcf . [6] ) . Another source of stress inhomogeneity inside the inclusions is a continuous variation of their 
mechanical properties such as in either cylindrically or spherically anisotropic particles (see, e.g., [48]). 
However, probably the most often investigated reason of such a stress inhomogeneity is an imperfect in- 
terphase (including sliding, debonding, cohesive phenomena, see for references, e.g., [6]). These interphase 
may represent weak intcrfacial layer due to imperfect bonding between the two phases and inter-diffusion 
and/or chemical interaction zones (with properties varying through the thickness and/or along the sur- 
face) at the interphase between the two phases. The thickness of interphase investigated usually ranges 
from h/a c = 0.01 for the conventional composites to h/a c = 2 for nanocomposites. The significance of 
interphase effects becomes important in nanocomposites due to their high surface-to-volume ratios. An 
alternative approach taking into account interfacial effect is based on the concept of surface stress and 
surface tension (see, e.g., [49], [50]). To the author's knowledge, in tens of publications dedicated to the 
influence of interphase on effective properties, the methods usually based on the hypothesis HI (such as, 
e.g., the Mori-Tanaka scheme and MEF) are exploited. Now all these estimations can be improved in the 
framework of the new approach as we did it in Figs. 2-10 [compare the results obtained for the initial 
B* [0] (x) and 12th B* [12] (x) iterations]. 

Other possible directions of successful applications of the proposed approach are three classes of 
problems where inhomogencitics of stress distributions in the inclusions arc generated by the nonlocal 
effects even for homogeneous ellipsoidal inclusions. The first two classes of these problems are described 
by both the special features of applied loading (statistically homogeneous media subjected to inhomo- 
geneous boundary conditions) and the special features of microstructure (FGMs, clustered materials, 
bounded media, contact of microinhomogencous media, macro-heterogencity insde microinhomogeneous 
medium, see for details and references [6]). In both cases, the known methods are based on the general 
integral equation (2.2) for the statistically inhomogeneous media when (r/)(y) ^const. (see for details 
and references [6]). However, Eq. (2.2) is just an approximation obtained from the exact Eq. (4.1) at the 
assumption (3.231). Using of more general Eq. (4.1) instead of the approximative Eq. (2.2) opens up great 
opportunities for detection of new effects in nonlocal micromechanics. The mentioned problems imply an 
estimation of nonlocal effective properties for composites through their constituents exhibit local consti- 
tutive properties. A new inverse problem is initiated by investigation of nanocomposites and formulated 
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as an estimation of local effective properties through the nonlocal mechanical properties of constituents. 
This problem was solved by Buryachenko [6] (see Section 18.2) in the framework of hypothesis HI. How- 
ever, it is well known in the context of micropolar elasticity that the strains are non-uniform even for the 
homogeneous elastic properties of the ellipsoidal inclusion subjected to the homogeneous remote loading. 
This sort of inhomogeneity is an encouragement for generalization of Eq. (4.1) to the composites which 
constituents are described by the nonlocal constituent laws. A subsequent step is the adoption of the 
new approach proposed in this paper for analysis of the generalized Eq. (4.1). However, more detailed 
consideration of nonlocal effects mentioned is beyond the scope of the current study and will be analyzed 
in other publications. 
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